% BL_OutpuGap_REStat_3 - Figures 6
%
% Replication files for the paper "Measuring the output gap with large datasets"
% by Matteo Barigozzi and Matteo Luciani, forthcoming on the Review of economics and statistics.
% 

% Written by Matteo Luciani, matteoluciani@yahoo.it

clear all; close all; clc; 
ML_graph_options
stampa=0; if stampa==1; disp('Results will be printed'); end

load RestatSetting                                                          % loads benchmark settings
v2struct(Setting);                                                          % "unpack" benchmark settings
[Y, Label, Name, Dates, cd, CBO] = ML_ReadDataHaver(Filename,trans,out);     % Load data in levels


%%% =============================== %%%
%%%  Initialize the model with PCA  %%%
%%% =============================== %%%
y=ML_diff(Y);                                                               % Data in 1st Differences
[T,N]=size(y);                                                              % size of the panel
[yy, my, sy]=ML_Standardize(y);                                             % Standardize
[f,lambda]=ML_efactors2(yy,q,2);                                            % estimate factor loadings, aka DFM on \Delta y_t
TT=(1:T+1)';                                                                % time trend
X=NaN(T+1,N); bt=X; b=zeros(N,2);                                           % preallocates variables for detrending
J=BL_TestLinearTrend(y); J(TR1)=1;                                          % Identify variables to be detrended
[X(:,J),bt(:,J),b(J,:)]=ML_detrend(Y(:,J));                                 % Detrend variables to be detrended
X(:,~J)=Y(:,~J)-repmat(mean(Y(:,~J)),T+1,1);                                % Demean variables not to be detrended
bt(:,~J)=repmat(mean(Y(:,~J)),T+1,1);                                       % ---------------------
b(~J,1)=mean(Y(:,~J))';b(~J,2)=0;                                           % ---------------------
if GDO==1 % --------------------------------------------------------------- % Restrictions for GDO
    lambda(1:2,:)=repmat(mean(lambda(1:2,:)),2,1);                          % same loadings
    b(1:2,2)=mean(b(1:2,2));                                                % same slope, different constant
    bt(:,1:2)=repmat(b(1:2,1)',T+1,1)+TT*b(1:2,2)';                         % same linear trend
    X(:,1:2)=(Y(:,1:2)-bt(:,1:2));                                          % Detrended GDP GDI
    sy(1:2)=repmat(mean(sy(1:2)),1,2);                                      % same std for standardization
end        % -------------------------------------------------------------- %
Z=X./repmat(sy,size(X,1),1);                                                % Standardize detrended variables  
F=Z*lambda/N;                                                               % Factors in levels as in BLL
if strcmp('VAR',model); [A0,v,AL]=ML_VAR(F,p,1);                            % Estimate VAR
elseif strcmp('VECM',model); [A0,~,~,~,~,v,~,AL]=BLL_VECM(F,p,d,10,1);     	% Estimate VECM
end
xi=Z-F*lambda';                                                             % idiosyncratic component

%%% ============================ %%%
%%%  Estimate the model with EM  %%%
%%% ============================ %%%
Z2=(Y-repmat(b(:,1)',size(Y,1),1))./repmat(sy,size(Y,1),1);                 % BL standardizations: (Y_t - a), i.e. approx. centered levels, divided by sy 
b1=b./repmat(sy',1,2); b1(:,1)=0;                                           % divide slope by sy, constant=0;
[A,Q,Lambda,R,F00,P0s,mu,type2]=ML_DynamicFactorSS_GDO_TV...                % State-space representation
    (AL,lambda(:,1:q),v,q,s,xi,F,rr,I0,I1,b1,A0,TV,star);                   % --------------------------
[xitT,PtT,PtTm,~,~,Ptt,~,A,Lambda,R,Q,b2]=ML_DynamicFactorEM_GDO_TV...      % EM-Algorithm
    (Z2,F00,P0s,A,Lambda,R,Q,s,q,d,p,mu,type2,model,maxiter,tresh,cc,GDO);  % ------------
T2=size(Z2,1); start=1;
BT=TT*b2';                                                                  % ML estimates of linera trend
isTV=find(type2==10)+max(p*q,q*(s+1));                                      % identifies TV coefficients
for ii=1:length(isTV) % --------------------------------------------------- % TV slopes or means
    id=find(Lambda(:,isTV(ii)));
    BT(:,id)=xitT(:,isTV(ii))*ones(1,length(id));                     
end                   % --------------------------------------------------- %
FF=xitT(:,1:q); FF1=[]; for ss=0:s;FF1=cat(2,FF1,FF(s-ss+1:end-ss,:));end   % Store ML Factors
for ss=1:s+1; lambda2(:,:,ss)=Lambda(:,(ss-1)*q+1:ss*q); end                % Store ML loadings
L=reshape(lambda2,N,q*(s+1));

FF2=[]; for pp=1:p;FF2=cat(2,FF2,FF(p-pp+1:end-pp,:));end                   % compute common shocks
rawshock=FF(p+1:end,:)-(A(1:q,1:q*p)*FF2')';                                % ---------------------

%%% =================== %%%
%%%  Common Components  %%%
%%% =================== %%%
T3=length(FF1(cc:end,:));
SY=repmat(sy,T3,1);
MY=repmat(b(:,1)',T3,1);
start2=start+s+cc-1;
Y2=Y(start2:end,:);
chi=(FF1(cc:end,:)*L'+BT(start2:end,:)).*SY+MY;                             % common component
zeta=Y2-chi;                                                                % idiosyncratic component    

%%% ================================= %%%
%%%  Estimate trend with PCA on LRCV  %%%
%%% ================================= %%%
[V, D]= eig(cov(FF(cc:end,:)));                                             % eigenvalue and eigenvenctors of VCV
[~, t2]=sort(diag(D),'descend'); V=V(:,t2); D=D(t2,t2);                     % sort eigenvalues and eigenvectors
Vt=V(:,1:q-d); Vc=V(:,q-d+1:q);
VVt=Vt*Vt'; VVc=Vc*Vc';


for jj=1:2
    if jj==1        % ----------------------------------------------------- % PCA
        FFt = FF(cc:end,:)*VVt;                                             % TREND Factors
        FFc = FF(cc:end,:)*VVc;                                             % CYCLE Factors     
    elseif jj==2    % ----------------------------------------------------- % RW with filter
        FF2=FF*Vt;                                                          
        SS=BL_StateSpaceLLLM(FF2,var(rawshock*Vt));        
        LLLM=BL_LLLM(FF2,SS,0,.1,cc,1);      
        FFt=LLLM.xitT(cc:end,1) * Vt';
        FFc=FF(cc:end,:)-FFt; 
        tau=[FF2 LLLM.xitT(:,1)];
    end             % ----------------------------------------------------- %
    FFt1=[]; FFc1=[];
    for ss=1:s+1
        FFt1=cat(2,FFt1,FFt(s-ss+2:end-ss+1,:));                            % lagged trend factors
        FFc1=cat(2,FFc1,FFc(s-ss+2:end-ss+1,:));                            % lagged cycle factors
    end
    chit(:,:,jj) = (FFt1*L'+BT(start2:end,:)).*SY+MY;                       % common stochastic trend plus linear trend
    chic(:,:,jj) = (FFc1*L').*SY;                                           % common temporary
end




%%% ============= %%%
%%% CBO Estimates %%%
%%% ============= %%%
J=ismember(CBO(:,1),Dates(start2:end));
cbot=100*log(CBO(J,2));                                                     % Potential Output CBO
cboc=Y2(:,1)-cbot;                                                           % Output gap CBO
cbonr=CBO(J,3);                                                             % Natural Rate CBO
cboug=Y2(:,75)-cbonr;                                                        % Unemployment gap CBO


%%% ========================= %%%
%%%  Loads benchmark results  %%%
%%% ========================= %%%

load RestatBenchmark2;
alpha1=RESTAT.alpha(1);
alpha2=RESTAT.alpha(2);

ogap(:,1)=RESTAT.chic(:,1);
pot(:,1)=RESTAT.chit(:,1);
ogap(:,2)=chic(:,1,2);
pot(:,2)=chit(:,1,2);


%%%%%  ======================================  %%%%%
%%%%%  ====  --------------------------  ====  %%%%%
%%%%%  ====  ---                    ---  ====  %%%%%
%%%%%  ====  ---  Plotting Results  ---  ====  %%%%%
%%%%%  ====  ---                    ---  ====  %%%%%
%%%%%  ====  --------------------------  ====  %%%%%
%%%%%  ======================================  %%%%%

Dates2=Dates(start2:end); Dates2d=Dates2(2:end); Dates4q=Dates2(5:end);     % New Dates for plots
j0=find(year(Dates2)==1970,1);                                              % Strarting point for all graphs 
nomefile='Restat_00bw_';
lg={'CBO','BL','BL$_{rw}$','$\widetilde{\mathrm{BL}}_{rw}$',...             % labels for legend
    '$\widetilde{\mathrm{BL}}_{rw}^0$'};                                    % -----------------
  
LS={'k-','k--','k:','k-.'};
DD=Dates2(j0:end);


%%% ========== %%%
%%%  Figure 6  %%%
%%% ========== %%%

ZZ=[cboc ogap];   ZZ=ZZ(j0:end,:);                                          % Output Gap - levels
ZZb{1}=BL_Band(ZZ(:,2),(RESTAT.ogapB(j0:end,:)),alpha1);                    % -------------------
ZZb{2}=BL_Band(ZZ(:,2),(RESTAT.ogapB(j0:end,:)),alpha2);                    % -------------------
SizeXY=ML_SizeXY(DD,ML_minmax([ZZ ZZb{1}]),1,5);                            % -------------------
SizeXY(1:2)=SizeXY(1:2)+[-60 60];                                           % -------------------
BL_ShadowPlotBW(ZZb,ZZ,DD,lg,SizeXY,LS([2 1 3]))
if stampa==1; print('-depsc','-painters','-r600',[nomefile 'OG_rw']); end      % -------------------
title('Output Gap - Level','fontweight','bold','fontsize',16)               % -------------------

ZZ=[cbot pot]; ZZ=ML_diff(ZZ(j0:end,:),1);                                  % Potential Output - 4Q
ZZb{1}=BL_Band(ZZ(:,2),(ML_diff(RESTAT.potB(j0:end,:),1)),alpha1);          % ---------------------
ZZb{2}=BL_Band(ZZ(:,2),(ML_diff(RESTAT.potB(j0:end,:),1)),alpha2);          % ---------------------
SizeXY=ML_SizeXY(Dates4q(j0:end),ML_minmax([ZZ ZZb{1}]),1,5);               % ---------------------
SizeXY(1:2)=SizeXY(1:2)+[-60 60];                                           % ---------------------
BL_ShadowPlotBW(ZZb,ZZ,Dates4q(j0:end),lg,SizeXY,LS([2 1 3]))
if stampa==1; print('-depsc','-painters','-r600',[nomefile 'PO_rw']); end      % ---------------------
title('Potential Output - 4Q% changes','fontweight','bold','fontsize',16)   % ---------------------  
    


 